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Abstract 


A  technique  for  determining  the  orbital  element  set  of  a 
sunlight-illuminated  object  detected  by  an  overhead  platform 
(when  passing  through  the  sensor's  field  of  view)  is  devel¬ 
oped.  The  technique  uses  a  Gauss  orbit  determination  tech¬ 
nique  to  find  an  initial  target  state  estimate  and  then  the 
estimate  is  refined  via  a  batch  weighted  least  squares  estima¬ 
tion  routine.  A  six  element  state  vector  consisting  of  three 
position  and  three  velocity  components  describe  the  state  at 
epoch.  It  was  found  that  the  Gaussian  method  produced  rea¬ 
sonable  initial  orbits  when  the  data  bias  was  sufficiently 
zero.  Each  analyst-supplied  slant  range  fit  the  data  equally 
well,  indicating  that  orbit  determination  is  impossible  with 
a  sinyla  set  of  data.  A  unique  series  of  events  where  the 
same  object  was  tracked  four  consecutive  days  was  fit  using 
the  developed  algorithm,  producing  favorable  results.  The  re¬ 
sults  of  two  single  data  set  events  and  one  multiple  collec¬ 
tion  events  are  presented. 
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Orbit  Determination  of  Sunlight-Illuminated  Objects 
Detected  by  Overhead  Platforms 


Due  to  the  multitude  of  objects  in  the  geostationary 
belt,  overhead  platforms  are  being  saturated  by  reflected 
sunlight  from  orbiting  objects  passing  through  the  sensors' 
field  of  view.  These  objects,  known  as  fastwalkers,  are 
creating  a  suspicion  that  some  uncatalogued  objects  may  exist 
or  are  being  cross-tagged  within  the  data  base.  The  North 
American  Aerospace  Defense  Command  (NORAD)  tasked  the  Foreign 
Technology  Division,  Flight  Performance  Division  (FTD/SQDF) 
to  analyze  these  15  to  30  minute  data  tracks  and  determine  the 
element  set,  identifying  the  object. 

The  objective  of  this  thesis  is  twofold;  to  determine  the 
feasibility  of  determining  the  element  set  of  an  orbiting 
object  from  space-based  metric  data  and  if  so,  perform  a 
commonality/occurrence  frequency  study  of  a  year's  worth  of 
collected  data.  The  project  can  be  expanded  to  warn  the 
sensor  operator  when  the  next  occurrence  will  be  so  that  he 
can  take  preventative  measures  to  protect  the  sensor  from 
damaging  cell  over-saturation. 


In  order  to  start  the  analysis,  a  few  underlying  assump¬ 
tions  must  be  made.  First,  the  sensor  ephemerides  are  known 
exactly,  since  this  is  the  best  true  baseline  information 
available.  Also,  sensor  location  is  a  type  of  Q-parameter, 
where  Q-parameters  are  defined  by  Day  as  those  parameters 
which  effect  the  observations  but,  for  some  reason,  cannot  be 
estimated  (5:3-1).  Other  examples  of  Q-parameters  are  atmo¬ 
spheric  density,  data  biases,  and  the  like.  P-parameters ,  on 
the  other  hand,  are  those  parameters  which  can  be  estimated 
from  the  given  data. 

Secondly,  the  objects  (targets)  are  assumed  to  be  non¬ 
thrusting  bodies  within  1000  km  from  the  sensor  since  the  col¬ 
lected  data  tracks  can  be  as  long  as  30  minutes,  implying  that 
the  object's  orbital  speed  is  nearly  the  same  as  the  sensor's 
orbital  speed.  This  assumption  along  with  the  fact  that 
orbital  perturbations  will  have  no  visible  effect  on  the 
target  during  the  span  of  the  data  track,  the  two  body  equa¬ 
tions  of  motion  will  be  sufficient  to  estimate  the  orbital 
element  set. 

The  simulated  data  for  this  project  were  generated  by  the 
Modularized  Vehicle  Simulation/Trajectory  Reconstruction 
Program  (MVS/TRP),  a  batch  weighted  least  squares  estimator 
originally  developed  to  validate  the  guidance  equations  of 
various  space  boosters.  The  data  were  generated  with  l-a 
Gaussian  white  noise  at  one  point  every  ten  seconds. 
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Table  I  Description  of  Simulated  Test  Case  Data 


Data  Rate:  lpt/10  sec 
Nominal  l-o  Data  Errors 
Azimuth:  3.5  X  10*’  Elevation:  4.3  X  10 


Research 

Since  the  fastwalker  data  base  goes  back  as  far  as  1972, 
it  is  safe  to  assume  that  the  problem  has  been  in  existence 
since  then.  It  is  most  prevalent  with  a  certain  ballistic 
missile  early  warning  satellite  system  located  at  various 
geostationary  locations.  Wong  in  his  paper  performed  his 
analysis  using  various  intensity  models  and  a  least  squares 
estimator  (MVS/TRP)  to  find  that  there  is  a  family  of  least 
squares  solutions  in  r  which  satisfy  the  collected  azimuth  and 
elevation  data  (21:23-25).  Little  was  done  to  resolve  this 
family  into  a  single  possible  fastwalker,  hence  his  conclusion 
was  that  the  problem  was  unsolvable. 

Some  undocumented  simulations  were  performed  in  the  Space 
Surveillance  Center  within  the  NORAD  Cheyenne  Mountain  Complex 
where  the  sensor  was  boosted  into  a  higher  orbit  and  then 
circularized  500  km  above  its  original  position.  Then  based 
on  the  viewing  angle  a  collection  of  other  satellites  in  its 


field  of  view  were  listed  at  a  particular  epoch  time.  This 
proved  to  be  futile  in  singling  out  one  particular  satellite 
as  the  fastwalker.  Other  than  the  above  mentioned  instances, 
no  other  work  has  been  performed  on  fastwalker  analysis. 

Approach 

Since  pure  angular  data  is  given  along  with  the  sensor 
ephemeris  over  the  data  span,  the  Gauss  orbit  determination 
method  will  be  used  to  determine  an  intial  orbit.  The  initial 
target  (fastwalker)  state  vector  can  then  be  refined  with  a 
batch  weighted  least  squares  estimator  to  further  update  the 
orbit. 

The  ephemeris  is  reported  in  latitude,  longitude,  alti¬ 
tude  and  time  since  no  velocity  information  is  required  to  fit 
the  data  (i.e.  the  data  is  also  position  related.  No  doppler 
or  other  velocity  related  effects  need  to  be  taken  into  ac¬ 
count  to  correct  the  data  model).  However,  an  accurate  veloc¬ 
ity  component  of  the  sensor  state  vector  roust  be  determined 
if  a  six  element  state  vector  will  be  propagated  over  a  few 
days. 

Assuming  the  first  ephemeris  position  data  point  is 
exact,  the  remainder  of  the  position  ephemeris  can  be  fit 
using  a  least  squares  algorithm.  The  first  difference  is  an 
excellent  first  guess  to  the  velocity  required  to  "hit”  all 
of  the  remaining  position  ephemeris  points.  Once  converged. 
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the  sensor  state  vector  will  be  the  best  possible  estimate 
given  the  data  at  hand. 

The  next  step  is  to  refine  the  derived  state  vector  from 
the  Gauss  algorithm  using  a  least  squares  estimator.  Assuming 
the  Gauss-derived  estimate  is  reasonable,  the  target  vector 
can  be  Improved  up  to  the  error  of  the  data.  The  1-a  error 
figures  of  merit  from  the  sensor  specifications  are  used  as 
an  input  into  the  program.  These  inputs  form  the  data 
weighting  matrix  required  by  the  least  squares  algorithm  to 
normalize  the  data  entries  (see  p.  37). 

With  the  fastwalker  orbit  now  defined,  it  can  be  propa¬ 
gated  forward  to  determine  the  next  encounter  with  another 
sensor  in  this  or  another  constellation.  This  is  easily  done 
with  the  Air  Force  Jet  Propulsion  Laboratory's  Long  Term  orbit 
Predictor  (LOP).  The  accuracy  of  this  predictor  has  been 
proven  with  simulations  matching  probe  flight  and  highly 
eccentric  orbit  trajectories. 
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II. 


This  chapter  presents  the  equations  necessary  to  accu¬ 
rately  model  the  motion  of  an  orbiting  object  over  short  and 
long  periods  of  tine.  This  includes  the  differential  equa¬ 
tions  of  motion,  transformation  from  an  inertial  state  vector 
to  classical  orbital  elements.  Aj.so  included  are  the  trans¬ 
formation  from  the  sensor  coordinate  system  to  the  inertial 
coordinate  frame,  which  is  the  computational  coordinate  sys¬ 
tem.  These  transformations  will  be  necessary  for  the  weighted 
least  squares  (also  known  as  differential  correction)  process 
where  the  observed  measurements  are  compared  to  the  assumed 
or  estimated  measurements. 

Coordinate  Systems  and  Transformations 

Many  different  coordinate  systems  are  quite  useful  for 
expressing  the  orbit  of  an  artificial  satellite.  Two  major 
coordinate  systems,  the  sensor  coordinate  frame  and  the  iner¬ 
tial  coordinate  frame  will  be  described  here. 

The  computational  coordinate  system  for  this  algorithm, 
the  frame  in  which  the  equations  of  motion  are  expressed  and 
integrated,  is  the  Earth  Centered  inertial,  or  ECI  coordinate 
system  shown  in  Figure  1,  p.  7.  The  origin  of  this  frame  ic 
at  the  center  of  the  ellipsoidal  earth  model;  the  X-Y  plane 
is  the  equatorial  plane  of  the  earth  and  the  Z-axls  is  coinci- 
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dent  with  the  earth's  axis  of  rotation.  This  coordinate 
system  is  fixed  in  inertial  space  with  X  pointing  toward  the 
first  point  of  Ares  at  some  arbitrarily  chosen  time  t^.  For 
this  program^  t,  is  the  midpoint  of  the  data  set  for  the  ini¬ 
tial  estimate  of  the  state  vector  and  the  time  of  the  first 
data  point  for  the  WLS  estimation  process. 

Two  latitude  types  are  normally  associated  with  an  ellip¬ 
soidal  earth  model  —  the  geocentric  and  geodetic  latitudes 
(see  Figure  2,  p.  9).  The  geocentric  latitude  is  defined  as 
the  angle  between  the  equatorial  plane  and  the  vector  from 
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the  center  of  the  earth  to  the  satellite  subpoint.  Geodetic 
latitude  on  the  other  hand  Is  the  angle  between  the  equatorial 
plane  and  the  line  perpendicular  to  the  subpoint  tangent 
(local  horizontal).  The  two  are  related  by  (1:97) 

tan  =  (a/b)*  tan  <p  (2.1) 

where  «p  Is  the  geocentric  latitude,  (p*  Is  the  geod-Jtlc  lati¬ 
tude  and  a  and  b  are  the  earth's  semlmajor  and  semlnlnor  axes 
respectively.  Geocentric  latitude  Is  used  for  all  computa¬ 
tions.  If  the  latitudes  are  Input  as  geodetic,  they  are 
converted  to  geocentric. 

The  other  coordinate  system  Is  the  sensor-fixed  system 
where  the  x-axls  is  pointing  toward  the  origin  of  the  ECI 
frame;  the  y-axls  Is  parallel  with  the  equatorial  plane  of 
the  earth  and  the  z-axls  forms  a  right-handed  coordinate 
system  where  x  X  y  =  z  (upper  case  letters  denote  the  Iner¬ 
tial  frame,  lower  case  letters  denote  sensor  frame). 

Within  the  sensor  frame,  the  target  Is  measured  In  azi¬ 
muth  and  elevation.  Azimuth  is  measured  In  the  detector  plane 
clockwise  from  "north”;  elevation  is  measured  up  from  the 
nadir  (x-axis)  direction  (see  Figure  3,  p.  11).  Both  can  be 
expressed  In  the  ECI  frame 

R.  Ry  P,  Ri  Pi 

cos  (El)  =  -  (2.2) 

R  P 
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6356.782  K" 


Eccentricity:  0,08182 


Figure  2 


Geoceuiric  Latitude 
Geodetic  Latitude 


Ellipsoidal  Earth  Model 


tan  (Az)  = 


sin  A 


cos  A 


sin  A  =  R  (R.  p,  -  Ry  p.) 
cos  A  =  pZ  (R/  +  R/)  -  R,  (R,  p,  +  Ry  Py) 


(2.3) 


(2.4) 

(2.5) 


where  R.,  Ry  and  R,  are  the  inertial  sensor  vector  elements  and 
p^,  Py  and  p,  are  the  inertial  slant  range  vector  elements. 
Projecting  the  range  unit  vector  upon  the  sensor  coordinate 
frame  produces  the  direction  cosines  L„,  Ly,  and  L„  where 


L„  =  cos  az 


(2.6) 


Figure  3 


Fastwalker  Sensor  Geometry 


Ly,  =  -sin  el  sin  az  (2.7) 

L,.  =  sin  el  cos  az.  (2.8) 

A  bias  and  scale  factor  is  added  to  each  azimuth  and  elevation 
calculation  to  account  for  other  metric  data  formats. 


The  two  coordinate  systems  are  related  by  two  rotations. 
The  first  rotation  is  about  the  sensor  y-axis  by  the  latitude 
and  then  about  the  sensor  z-axis  by  the  local  hour  angle. 
Mathematically, 


p.l 

COS  a  cos 

B 

sin  a 

-cos  a  sin 

B 

Pn 

Pyi 

S 

-sin  a  cos 

B 

cos  a 

sin  a  sin 

B 

Pf. 

Pti 

sin  B 

0 

cos  B 

P» 

where 

p.i,  etc.  *  inertial  frame  slant  range  vector  elements 
p„,  etc.  *  sensor  frame  slant  range  vector  elements 
a  *s  -local  hour  angle  ■  -[  +  w,(t+to)  +  a,  ] 

«  Greenwich  hour  angle  at  midnight 
w,  «=  earth  ^s  rotation  rate 
t+to  “  time  since  midnight 
=  sensor  longitude 
B  s  -V,  »  sensor  latitude. 

Geopotential  and  Gravity  Determination 

Assuming  a  non-thrusting  satellite  above  the  earth's 
atmosphere,  the  only  significant  force  acting  on  the  body  over 
a  data  span  is  gravity.  Lillard  in  his  thesis  assumes  negli¬ 
gible  longitude  dependence,  hence  the  earth's  geopotential 
becomes  (16;i8) 
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U  =  n/r  [  1  -  S  J„  (a./r)"  P„(sin(«p))  ], 


(2.10) 


^  =  Earth's  gravitation  constant  398600.8  kroVsec’ 
r  =  distance  from  center  of  the  earth  to  the  satellite 
J„  =  harmonic  constants 

a.  =  earth's  equatorial  radius,  6378.135  km 
P„  =  associated  Legendre  function 
(p  =  geocentric  latitude. 

This  is  also  known  as  the  zonal  potential,  which  will  be  quite 
sufficient  for  the  initial  orbit  estimate  due  to  short  data 
tracks  relative  to  the  time  for  perturbations  to  take  effect. 
Hence  four  expansion  terms  are  adequate.  For  long  term  propa¬ 
gation,  tesseral  and  sectoral  harmonics  and  other  effects  such 
as  third  body  effects,  solar  drag,  etc.  must  be  included. 

The  acceleration  due  to  gravity  is  derived  from  the  above 
potential.  Since  U  is  a  solution  to  Laplace's  equation,  the 
gravity  g„  is  simply  the  gradient  of  U,  or 

g„  =  -  V  (XJ).  (2.11) 


By  setting 


sin(<p)  =  2/r 

the  geopotential  and  gravitational  acceleration  can  be  written 
entirely  in  the  computational  coordinate  system,  eliminating 
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the  need  to  transform  coordinates.  Therefore,  the  gravita¬ 
tional  acceleration  becomes 


where  ('•)  denotes  the  second  derivative  with  respect  to  time. 
These  equations  of  motion  are  integrated  using  a  Runge-Kutta- 
Nystrom  numerical  integrator  to  propagate  the  state  along  in 
time . 
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For  this  specific  problem,  Escobal  states  the  Gaussian 
method  of  orbit  determination  is  "second  to  none"  when  the 
time  between  the  observations  is  small  (6:272).  The  Laplace 
method  was  found  to  be  unsuitable  since  the  direction  cosine 
rates  and  accelerations  were  small  enough  to  produce  near¬ 
singular  matrices. 

The  Gauss  method  hinges  on  the  fact  that  only  two  linear¬ 
ly  independent  vectors  are  needed  to  define  the  orbit  plane, 
assuming  negligible  perturbation  effects  during  the  data  span. 
Therefore  a  third  vector  can  be  written  as  a  linear  combina¬ 
tion  of  the  independent  vectors.  Thus  it  is  possible  to 
determine  a  set  of  constants  a,  b,  and  c  not  all  zero  such 
that 

ari  +  br,  +  cr,  *  0.  (2.15) 

Arbitrarily  choosing  r,  as  the  dependent  vector  and  redefining 
constants  to  c,  and  c, 

r,  =  c,r,  +  c,r,.  (2.16) 

Crossing  r,  with  r,  and  r,  three  parallel  vectors  perpen¬ 
dicular  to  the  orbit  plane  are  found  whereby 

r,  X  r,  =  r,r,sin  »/„  *  2A„lf 

r,  X  r,  =  r,r,sin  =  2A„W  (2.17) 

r,  X  r,  *  rir,sin  i/j,  *  2Ai,W 

14 


where  are  the  areas  of  the  triangles  formed  by  the  respec¬ 
tive  radius  vectors  and  i/i,  is  the  angle  between  the  respective 


vectors . 
Substituting, 


AjjW  —  CjAijM  , 


CjAj,ff 


or 


(2.18) 


The  Gauss  method  only  requires  observations  at  three 
different  times.  Since  data  sets  of  more  than  three  observa¬ 
tions  are  available,  the  set  to  be  analyzed  roust  be  "averaged" 
in  some  sense  (in  this  case,  time).  The  number  of  data  points 
was  divided  into  three  batches  and  the  average  time  of  each 
batch  was  used  as  the  batch  tiroes.  An  n'^-order  least  squares 
polynoroial  (determined  by  the  using  analyst)  is  fit  to  each 
data  coordinate  to  smooth  through  the  noise.  Then  based  upon 
the  resulting  equation,  ECl  direction  cosines  and  sensor 
location  for  each  batch  time  can  be  calculated  and  stored. 

The  above  area  ratios  can  be  expressed  as  a  function  of 
time.  For  instance. 


A„  =  l/2lf  •  (r.  X  r,] 


(2.19) 
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along  with  expressed  in  terns  of  r,  in  £  and  g  series  (see 
p.  28  for  derivation),  i.e. 


r,  =  f»r,  +  g.r'. 


Substituting, 

Aj,  =  1/211  ((fir,  +  gir',)]  X  r,]  =  -hgi/2 

where  h  is  the  orbital  angular  momentum.  Similarly  A, 
expressed  as 

A„  =  1/2W  [(f,r,  +  g^r',)]  X  r,]  =  -hg,/2 

Note  that  there  a  common  area  A,,  where 

Ai,  =  1/2W  •  [(f,r,  +  g.r',)]  X  (f,r,  +  g,r',)] 

or 

Ai,  =  l/2h  (fig,-f,gi) 

In  terms  of  the  f  and  g  expansions,  the  coefficients  c 
are 


g. 

Cl  =  - 

fig,-f,g» 

g. 

c,  - - 

f.g,-f,gi 


(2.20) 

(2.21) 
,  can  be 

(2.22) 

(2.23) 

(2.24) 
ii  and  Cj 

(2.25) 
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Truncating  the  f  and  q  series  after  the  second  tern,  c,  and  c, 
can  be  approximated  by 


Cj  a 


C,  a  — 


r. 


1  +  -  (r,.'  -  r/) 

6 


u, 

1  +  -  (’‘i/  -  t/) 

6 


(2.26) 


where  u,  =  ^i/r^  By  allowing 


A,  = 


A,  =  - 


n 


B,  = 


B,  =  - 


6t,, 

Ti 


U 


6t„ 

Cl  a  Ai  +  BiU,  and  Cj  »  A,  +  B,u,,  where  c. 


(Ti/  -  T,') 

-  r,‘), 

-1. 


The  vector  equation  r,  =  Cir,  +  c,rj  can  be  modified  using 
the  basic  relation  r  =  p  -  R  to  produce  the  equation 

CiP,  +  c,p,  +  c,p,  =  c,R,  +  c,R,  +  c,R,. 

The  vector  pi  can  be  written  as  a  product  of  the  direction 
cosine  matrix  and  magnitude  p.  The  resulting  equation 
takes  the  form 


CiP,L,  +  c,p,L,  +  c,p,L,  =  G 


(2.27) 
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where  G  =  CiRi  +  c^,  +  CjRj.  Thus  the  previous  equation  can  be 


expressed  in  matrix  form  as 


byj  Lyj 

- p 

0  p 

_ J 

G. 

Gy 

(2.28) 

L., 

L.a  L., 

CjPj 

G. 

It  is  possible  to  form  the  inverse  of  the  matrix  L, 
solving  for  the  c^p^  matrix.  If  a^,  denotes  the  individual 
elements  of  L*'  and  vectors 


A  =  [A,,  -1,  A,]"  ,  B  =  [B,,  0,  B,]"  ,  X  =  [X,,  X,,  , 

V  =  [Y,,  Y,,  YJ"  ,  and  Z  =  [Z,,  Z,,  Z,y 


then  the  slant  range  Pt  for  each  batch  time  is 


Ai*  +  B/u, 

Pi  -  - 

Ai  +  BiU, 

p,  =  A,*  +  B/u’  (2.29) 

A,*  +  B/u, 

P)  =  - 

A,  +  B,u, 


where 


A,*  —  (SjiA'X  +  aijA'Y  +  a^jAZ) 
B/  =  (a^BX  +  a^jBY  +  a„BZ) 
A,'  =  -(a,iAX  +  a,jA  Y  +  a„A  Z) 
Bj*  =  -(ajiB  X  +  a„B  Y  +  a,,B  Z) 
A/  =  (a,, AX  +  a„AY  +  a„AZ) 
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B,*  =•  (a^BX  +  a„B‘Y  +  a„BZ). 

To  determine  the  magnitude  of  the  target's  radius  vector 
another  equation  is  required.  Another  linearly  independent 
relation  can  be  derived  by  dotting  the  vector  equation  relat> 
ing  target  position  to  sensor  position  into  itself,  producing 

r,’  =  P,'  -  2pL,R,  +  R,’.  (2.30) 

The  scalar  products  LR  are  known  from  the  collected  angular 
observations  and  sensor  location.  At  the  center  batch  time 
tj,  let 


C,  =  -2(X,L„  +  +  Z,L.,). 

Substituting  p,  and  c,  into  the  resultant  dot  product  equation 

r/  -  (C,A/  +  A/*  +  R/)r/  -  u(C,B/  -*•  2A/B/)r,>) 

-  u’B/'  =  0.  (2,31) 

One  is  not  assured  of  a  real  root  since  this  equation  is 
of  even  order.  By  factoring  this  equation,  all  possible 
combinations  of  reasonable  roots  can  be  explored  by  the  ana¬ 
lyst.  The  classical  Newton-Raphson  root  finding  technique 
finds  only  one  root  based  upon  an  initial  guess.  Hence  the 
former  method  was  chosen. 

Velocity  is  found  by  first  taking  a  Taylor  series  expan¬ 
sion  of  r,,  i.e. , 
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eo 


Where  t,,  =  (t,  -  being  the  epoch  time.  When  the  expan¬ 

sion  is  written  for  batch  tines  1  and  3,  tine  2  being  epoch 


r,"  r,'"  r’'' 

r.  -  r,  =  -T„r/  +  r,/  -  -  t,/  -  +  t„‘  —  +  0(Ti,*)  (2.33) 

2  6  24 


^  <1  ««•  l-iv 

r,  -  r,  =  T„r/  +  t„’  -  +  t„’  -  +  t„*  —  +  0(t„»)  (2.34) 

2  6  24 


Multiplying  eguation  (2.33)  by  t„,  equation  (2.34)  by  tj,  and 
adding  produces  an  equation  void  of  rj' .  Now  multiplying  the 
equation  (2.33)  by  and  equation  (2.34)  by  t,,’  and  adding 
produces  an  equation  void  of  r,”.  Differentiating  each  of 
these  new  expressions  yields 


T.,r.  -  r„r,”  +  1/2  r/''  +  0(r/) 

-r„’r,"  +  (r„’  -  T./)r,"  +  r./r,"  -  r„T„r.,  r,’"  + 

0(r,'"). 


(2.35) 


(2.36) 


Assuming  the  0(r*’')  terms  can  be  ignored,  one  can  solve  for  r/'' 
and  r,"  and  substituting  the  respective  derivatives  into  the 
equation  void  of  r,”  derived  above  generates 


20 


<li  ■  Gi  +  -  ,  i-1,2,3  (2.38) 

r.’ 

and 

V  -  -dir,  +  d^,  +  d,r,.  (2.39) 

At  this  point,  the  orbit  is  considered  determined  since 
both  r  and  v  are  known. 
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Orbital  Element  Sets 


Orbital  elements  are  sets  of  six  or  more  parameters  which 
describe  an  orbit,  whether  be  In  size,  shape,  orientation, 
energy,  momentum,  etc.  of  the  three  main  types  —  classical, 
equinoctial,  and  Delaunay  —  the  equinoctial  element  set 
provides  the  best  accuracy  since  there  are  no  singularities 
as  there  are  In  the  classical  and  Delaunay  sets,  except  at 
1=180*.  The  drawback,  however,  is  that  the  equinoctial  set  is 
rather  abstract  compared  to  the  classical  set. 

The  equinoctial  elements  in  terms  of  the  classical  ele¬ 
ment  set  for  a  direct  orbit  are: 

a  =  a  A=Mo  +  w-»-n 

h  =  e  sin  (w  +  n)  k  =  e  cos  (w  +  0)  (2.40) 

p  =  tan  (i/2)  sin  n  q  =  tan  (i/2)  cos  n. 

Deriving  these  elements  in  terms  of  the  state  vector  re¬ 
quires  a  little  more  knowledge  of  the  geometry  involved. 
Cefola  in  Figure  4  (p.  23)  shows  graphically  what  the  equinoc¬ 
tial  set  represents  (4:2).  The  orientation  of  the  coordinate 
system  is  defined  by  an  angle  equal  to  the  longitude  of  the 
ascending  node,  O,  and  the  unit  vector  w  perpendicular  to  the 
orbit  plane.  Unit  vector  f  is  in  the  orbit  plane  specified 
by  i/2  down  from  Ares  and  an  angle  n  down  from  the  node. 
Vector  g  forms  the  right-handed  triad  such  that  f  X  g  =  w. 
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Figure  4  Equinoctial  Coordinate  Frame 


An  orbit ^8  shape  is  defined  by  the  elements  h  and  k, 
which  are  the  components  of  the  eccentricity  vector  projected 
upon  the  f  and  g  vectors.  The  element  pinpoints  the  posi¬ 
tion  of  the  satellite  on  the  orbit  since  it  is  the  sum  of  the 
three  angles  measured  in  the  plane.  The  final  two  elements, 
p  and  q,  are  the  most  abstract.  In  short,  these  elements 
merely  characterize  the  orientation  of  the  orbit.  Both  are 
required  in  the  rotation  matrix  to  transform  the  inertial 
frame  to  the  equinoctial  frame  (4:3). 


For  the  direct  orbit  case,  the  unit  vectors  can  be  mathe¬ 


matically  expressed  as 


1  f  1  -  p»  +  qM 

f  =  .j  2pq  \ 

1  +  p’  +q*  [  -2p  J 

1  f  2pq  1 

g  .  J  1  +  P>  -  q‘  I 

1  +  p’  +q*  [  2q  J 

1  f  2p  ] 

w  =  \  -2q  \ 

1  +  p^  +q'  [  1  -  p’  -  q’  J 


(2.41) 


(2.42) 


(2.43) 


Transformation  from  an  ECI  State  Vector  to  Equinoctial  Ele- 

mfinta 

The  transformation  from  the  state  vector  to  the  equinoc¬ 
tial  elements  is  a  fairly  straight  forward  once  the  coordinate 
frame  is  understood.  An  orbit's  specific  energy  directly 
produces  the  semimajor  axis  since 


a  =■  - 


(2.44) 


The  eccentricity  vector  is  given  by 


• 

* 

1  u 

1 

> 

_ 1 

r  -  (r  v)  v  ■ 

(2.45) 
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Another  way  to  express  the  unit  normal  to  the  orbit  plane  Is 
through  the  angular  momentum,  r  X  v.  Hence 

r  X  V 

e  =  -  .  (2.46) 

|r  X  v| 

However,  this  expression  has  the  same  exact  meaning  as  the 
previous  w  equation  In  terms  of  p  and  q.  Therefore,  p  and  q 
are 


p  »  -  q  =  -  -  (2.47) 

1  +  e,  1  + 

and  the  equinoctial  elements  h  and  k  are  merely 

h-eg,  k-ef.  (2.48) 

The  remaining  element  to  derive  Is  the  mean  longitude, 
Xo.  First  compute  the  position  coordinates  X  and  Y  relative 
to  the  orbital  frame  by 


X  -  X  f,  Y  -  x  g. 


Then  compute 


(1  -  k’6)X  -  hkBY 

cos  F  -  k  -»•  - 

a(l  -  h’  -  k’)‘^’ 


(1  -  k’6)Y  -  hkBX 

sin  F  »  h  +  - 

a(l  -  h’  -  k*)*'* 


(2.49) 


(2.50) 
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where  the  variable  B  is 


1 

fl  a - 

1  +  (1  -  h'  - 

Then  the  mean  longitude  is  given  by 

A,  *  F  “  k  sin  F  +  h  cos  F.  (2.51) 

These  equations  are  mainly  used  to  compute  the  numerical 
partials  for  the  Jacobian  transforming  the  ECI  state  vector 
covariance  matrix  to  the  equinoctial  element  covariance  ma¬ 
trix. 
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III.  Propagation  Methods 


in  order  to  calculate  the  sensor  and  target  state  vectors 
at  tines  other  than  the  epoch,  the  state  must  be  propagated 
via  some  numerical  technique.  There  are  various  methods  to 
propagate  the  state  along  in  time  so  that  the  calculated  mea¬ 
surements  can  be  compared  to  the  actual  observations.  For 
this  thesis,  however,  three  main  methods  are  used,  each  having 
a  particular  purpose.  Two  methods  are  used  to  propagate  the 
state  —  the  4"‘-0rder  Runge-Kutta-Nystrom  numerical  integrator 
and  f  and  g  series  expansions.  A  state  transition  matrix  is 
used  to  propagate  the  covariance  matrix,  which  consists  of 
partial  derivatives  of  the  present  state  with  respect  to  the 
estimated  epoch  state.  The  covariance  matrix  will  be  examined 
in  the  next  chapter.  Each  propagation  method  will  be  briefly 
described. 

4th  order  Runqe-Kutta-Nvstrom  rRKN)  Numerical  Integrator 

The  RKN  integrator  is  a  closed  form  predictor-corrector 
numerical  integrator  specifically  suited  for  2,^  order  differ¬ 
ential  equations.  One  of  the  novelties  of  this  integrator  is 
that  the  position  and  velocity  at  time  t+at  are  direct  out¬ 
puts,  so  only  one  integration  pass  is  required  per  time  step. 

Battin  does  an  excellent  job  of  deriving  the  RKN  inte¬ 
grator  in  minute  detail  (2:577).  However,  Kreysig  describes 
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the  method  in  terms  of  a  scalar  function  which  is  easier  to 


follow  (15:1078). 

The  RKN  numerical  integrator  algorithm  is 
X  =  Xa  +  hy^  +  h*(ka  +  2ki)/6  +  0(h*) 

y  =  y«  +  h(kc  +  4ki  +  k,)/6  +  0(h*) 


(3.1) 


where 


kj  —  f  (  ta  /  Xa  ) 

=  f(ta  +  h/2,  Xa  +  hya/2,  h'ka/8) 
ka  =  f(ta  +  h,  Xa  +  hya ,  hTc,/2 ) 

This  algorithm  is  best  suited  for  accurate  propagation 
of  a  state  vector  in  an  estimation  routine  since  a  precise 
state  estimate  is  required  at  each  data  point.  In  addition, 
other  effects  can  be  included  in  the  function  f(t,x,y)  to 
exactly  model  the  physics  of  the  problem. 

F  and  G  Series 

The  f  and  g  series  expansions  are  series  solutions  of  the 
two-body  equation  of  motion.  In  other  words,  any  position  r 
can  bo  expressed  in  terms  of  the  epoch  position  rg  and  the 
time  r  from  the  epoch  to  the  time  desired.  Mathematically, 
any  function  can  be  expressed  as  an  infinite  series,  i.e. 
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(3.2) 


Which  is  the  form  of  a  Taylor  Series  about  zero  (also  known 
as  a  MacLaurin  Series.) 

Assuming  that  the  above  series  is  a  solution,  differen¬ 
tiate  the  two-body  equation  of  motion  three  times  (making  this 
expansion  of  fifth-order)  holding  the  epoch  time  t^  constant 
and  substituting  ro”  =  urj  (u=u/r’)  wherever  ro"  appears.  The 
resulting  equations  are 


ro‘“  =  -u'ro  -uro' 

ro"  =  (-u«  +  u’)ro  -  2u'ro'  (3.3) 

r*''  »  (-u“‘  +  4uu')r,  -  (3u''  -  u’)ro'. 

Substituting  these  expressions  into  the  Taylor  series  produces 
the  familiar  equation 


r  =  fr,  +  gro' , 


(3.4) 


Where 

f  ■  1  -  1/2  ut’  -  1/6  u't’  -  1/24  (u”  -  u’)f*  -  l/120(u"‘ 

-  4UU')’’’  (3.5a) 


and 

g  ■  r  -  1/6  ur’  -  1/12  u'r*  -  l/120(3u''  -  u’)t‘.  (3.5b) 
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In  this  form  the  f  and  g  expressions  are  useless  from  a 
practical  point  of  view  since  higher  order  derivatives  of  u 
are  not  known.  Escobal  states  Lagrange  proceeds  as  follows 


(6:109).  Introduce  the  new  variables  p  and  q  such  that 


where 


so, 


r»p 


d(r*) 

dr 


r'q 


d’(r’) 

dr’ 


(3.6) 


rr’  r  •  r' 
r’  r’ 


.  /  J  _  —3 


r'u 


q  3 


r’ 


3^i  dr 
r*  dr 


3n  d(r’) 


2r’  dr 


1  d’(r’)  1  d(r’)  dr 

2r’  dr’  r’  dr  dr 


(3.7) 


1  d’(r’)  dr  1  d’(r’) 

q'  =  -  - - + - 

r’  dr’  dr  r’  dr’ 


which  become  by  Introduction  of  the  previously  defined  p  and 

q 


u'  =  -3up 

p'  =  q  -  p’  (3.P) 

q'  =  -(up  +  2pq) . 
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Differentiating  u'. 


u«  =  -3(up'  +  u'p) , 


which  becomes 

U”  =  -3[U(q  -  2p’)  -  3up*]  (3.9) 

via  substituting  p  and  q.  This  process  can  be  repeated  to 
eliminate  n' ,  p^  and  q'  from  each  derivative.  Therefore,  the 
5“  order  f  and  g  series  expansions  are 

f  =*  1  -  1/2  ur’  +  1/2  upr’  +  1/24  (3uq  -  15up'  +  u’)r‘  + 

1/8  (7up’  -  3upq  -  u’p)T*  (3.10) 

g  =  T  -  1/6  ut’  +  1/4  upr*  +  1/120  (9uq  -  45up’  +  u’)r*. 

(3.11) 

These  expressions  are  used  to  obtain  a  quick  initial 
state  estimate  at  the  beginning  of  the  data  track.  Since  the 
series  is  truncated  at  0(t*),  this  method  is  definitely  not 
suited  for  long  term  propagation.  As  long  as  the  initial 
target  state  estimate  is  "good”,  the  least  squares  estimator 
can  properly  refine  it. 
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state  Transition  Matrix 


In  order  to  create  the  partial  derivative  matrix  required 
for  the  least  squares  process  /  th^2  partial  derivatives  of  Az^ 
and  Ell  must  be  taken  with  respect  to  the  state  vector  at  the 
epoch  time.  Partial  derivatives  of  the  i“  data  point  can  be 
related  to  the  epoch  state  vector  by  the  state  transition 
matrix,  ♦(ti*i,ti).  This  matrix  is  derived  in  many  control 
theory  texts,  but  its  use  for  our  purposes  Wiesel  derives  it 
in  terms  of  the  required  dynamical  equations  (18:26).  In 
short,  the  state  transition  matrix  demonstrates  how  the  state 
changes  in  time  but  in  matrix  format. 

The  state  transition  matrix  is  derived  from  the  solution 
of  the  matrix  differential  equation 

d 

—  ♦(t,t,)  =  A(t)  4Kt,tJ  (3.12) 

dt 

where  in  this  case  ♦  is  more  appropriately  called  the  equa¬ 
tions  of  variation  since  it  contains  the  effects  of  all  tra¬ 
jectories  near  the  reference  trajectory.  In  other  words,  it 
is  not  linearized.  If  one  assumes  that  the  changes  in  each 
element  is  negligible,  «  is  a  constant  matrix.  The  matrix 
also  takes  on  some  interesting  properties: 

♦  (t,t)  =  X 

♦  (t,,tj  =  ♦(t,,tj  ♦(ti,tj 
♦  (t,tj  =  ♦(t„t)-'. 
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(3.13a) 

(3.13b) 

(3.13c) 


The  partial  derivative  matrix  of  the  data  with  respect 
to  the  state  vector  at  t^  can  be  related  to  the  epoch  state 
vector  Xe  by  the  chain  rule  (equation  3.13b).  Let  be  the 
partial  derivative  comparing  the  changes  In  the  data  to  the 
state  elements  at  time  ti.  Then  Hi  Is  related  to  the  epoch 
vector  by 

Hi/„  =  Hi  ♦(ti,ti.i)  «-(ti.i,ti.,)  ...  *(ti,ti)  ♦(tiitj  (3.14) 

where  Hi,,  is  the  partial  derivative  matrix  of  the  data  at 
time  ti  with  respect  to  the  epoch  state  vector  x.  which  is 
required  for  use  in  the  least  squares  estimator. 
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IV.  Refining  the  Initial  State.  Estimate 


The  nain  method  of  refining  the  state  vector  derived  from 
the  Gauss  orbit  method  is  least  squares.  The  method  of  least 
squares,  developed  by  Gauss  for  his  astronomical  studies  in 
1795,  is  concerned  with  estimating  the  values  of  a  set  of 
parameters  of  a  measurement  model,  which  relates  the  measure¬ 
ments  to  the  parameters  to  be  estimated  (16:35).  These  esti¬ 
mated  measurements  are  then  compared  to  a  set  observations  to 
produce  residuals.  The  sum-squared  residuals  are  minimized. 
The  set  of  estimates  which  produces  this  minimum  is  said  to 
be  optimal  in  a  least-squares  sense  (i.e.  the  azimuth  and 
elevation  data  are  the  observations).  The  initial  state  esti¬ 
mate  produces  a  corresponding  set  of  observations  which  are 
subtracted  from  the  measured  observations  to  produce  residu¬ 
als. 

The  parameters  to  estimate  are  called  the  state  variables 
at  some  fixed  time  t,  known  as  the  epoch  time.  If  we  let  the 
vector  represent  the  set  of  measured  azimuth  and  elevation 
at  time  t^,  the  observation  model  may  be  represented  as 

z  “Gj(X,,to)  +  ,  i  "  1,2,. ..,n  (4.1) 


where 

n  a  total  number  of  measurement  times 

Xo  *  the  initial  estimate  of  the  state  vector  at  epoch 
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Gi  •  the  vector  of  functions  relating  the  state  vector  to 
the  observations 

Cl  -  the  error  vector  of  the  observations  at  tine  t^. 


Therefore,  the  observations  can  be  notatlonally  expressed  as 


1 

-  AZi 

+ 

Ca.  t 

El.*t 

-  El, 

+ 

Cti  1 

where 

Ca>i/  Cii  1  ”  the  1*^  azimuth  and  elevation  measurement  error 

AZtfei,  Elw  1  *  the  1**  measured  azimuth  and  elevation 

AZi  and  El|  "  the  true  azimuth  and  elevation  at  time  t^. 

The  minimum  state  vector  for  an  orbiting  object  Is  the 
six-element  vector 

X.-  (X.,  y„  z„  X/,  y,',  z/)  (4.3) 

assuming  a  high  altitude  non-perturbed  satellite.  Perturba¬ 
tions,  thrusting  and  other  effects  add  extra  elements.  Since 
one  of  the  aseumptions  is  the  object  is  non-thrusting  and  the 
data  track  Is  small  with  respect  the  time  It  takes  for  pertur¬ 
bations  to  take  effect,  these  elements  are  Ignored. 

If  the  observation  relations  were  linear,  one  could  take 
any  .o  of  them  and  solve  for  the  unknown  parameters,  where  m 
Is  the  number  of  unknown  state  elements  in  the  relations. 
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However,  the  observation  relations  are  non-linear  which  forces 
us  to  expand  the  relations  in  an  infinite  series,  truncate  and 
then  iterate.  This  is  also  known  as  linearization. 

A  perfect  sensor  would  measure  the  perfect  state  x.  pro¬ 
ducing  the  exact  set  of  measurements  z..  Then  the  estimated 
state  can  be  represented  as  x  »  x»  ^x.  Hence  the  error 
vector  e  is  given  by 

€  =  z  -  z,  =  G(x*  +  6x)  -  H(x„) 

dG 

as  —  Sx  (4.4) 

dX 

Let  Xo  be  the  estimate  of  x*.  Then  the  first-order  ap¬ 
proximation  of  the  error  by  expanding  the  observation  rela¬ 
tions  in  a  Taylor  series  is 

dAZ, 

€*.1  =  AZoti  -  AZt  -  -  AX  +  O(AX’)  (4.5) 

dXo 

aEii 

Ct,  ,  =  Eloti  -  El,  -  -  AX  +  O(AX’)  (4.6) 

dXo 

Where  ax  is  the  difference  between  the  initial  state  and  the 
state  estimate.  Then  the  partial  derivative  matrix  of  the 
observation  relation  G  to  the  state  elements  is 
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0 


0 


0 


dG 

dX 


dAZj 

dAZi 

dAz, 

ax. 

3y, 

a  2, 

aEl, 

dEl( 

dEl, 

ax. 

1 

o 

32, 

0 


(4.7) 


When  the  higher  order  terms  are  ignored,  the  Taylor  series 
expansion  produces  a  linear  set  of  equations  In  ax.  Now  we 
can  apply  the  concept  of  non-linear  least  squares. 

Let  e  be  the  2n  X  1  matrix  of  error  elements  azimuth  and 
elevation  for  each  observation  time 


H 


ti  I 

/  1  r 

,  /  •  • 

.  ,  € 

dGi  dG] 

aG) 

aG„ 

ax,  ax. 

• 

X 

1  ^ 

•  ~ 

ax. 

AZ,b  1 

-  AZi(X, 

-t,) 

El*. 

-  E1](X. 

-t,) 

AZ,5  1 

-  AZ,(X, 

-t,) 

El,., 

-  E1,(X, 

-t,) 

to  = 


(4.9) 


(4.10) 


Az^,  -  Az,(X,t„) 

[  EU.  -  El.(x,tJ  J 


With  the  above  matrix  definitions,  the  linearized  observation 
relations  can  be  rewritten  as 


€  =  to  -  MAX. 


(4.11) 
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There  are  two  types  of  residuals  —  azimuth  and  eleva¬ 
tion.  Elevation  is  normally  more  accurate  than  azimuth,  so 
the  least  squares  process  will  try  "fit"  azimuth  before  It 
considers  fitting  elevation.  Therefore  weights  must  be  ap¬ 
plied  to  the  data  to  normalize  the  residuals  (7:48).  This  way 
the  least  squares  process  will  consider  each  observation 
equally.  Also  weighting  gives  us  the  opportunity  to  Ignore 
unwanted  or  erroneous  data. 

Assuming  that  the  data  is  corrupted  by  Gaussian  white 
noise  (each  observation  Is  uncorrellated  In  time  and  space 
with  zero  mean),  the  weighting  matrix  W  is  a  diagonal  matrix 
with  the  weightings  on  the  diagonal.  This  is  not  true  since 
a  sensor  takes  raw  measurements  and  then  transforms  them  into 
something  which  can  be  understood.  Also,  If  the  observations 
are  measured  by  a  moving  sensor  they  are  correlated  in  time. 
However,  these  effects  are  generally  small  and  the  assumption 
Is  justified  (5:B-3). 

It  is  clear  when  weighting  is  introduced,  the  estimate 
will  be  a  direct  function  of  the  weights.  Lillard  states  the 
optimal  weighting  matrix  is  the  inverse  covariance  matrix  of 
the  measurement  errors  (16:43).  Considering  the  above  assump¬ 
tions,  the  weighting  matrix  is  composed  of  error  variances  on 
the  diagonal  with  the  variances  of  each  observation  data  type 
being  equal. 
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Introducing  the  2n  X  2n  weighting  matrix  the  matrix 


error  equation  is 


(4.12) 


Where  the  v  subscript  indicates  the  weighted  error  equation. 
"Squaring”  e  produces  the  squared  residual  matrix 

c/c,  -  -  w‘'"Hax)  (4.13) 


Now  taking  the  partial  derivatives  with  respect  to  /^x  and 
since  we  are  looking  for  a  minimum,  setting  it  equal  to  zero 
gives  the  standard  weighted  least  squares  difference  equation 


AX  -  (H’WM)**  H’Wto 
F  -  (H’WH)*'. 


(4.14) 


with  being  the  covariance  matrix  associated  with  the  esti¬ 
mate  where  the  diagonal  terms  are  the  variances  and  the  off- 
diagonal  terms  are  the  covariances. 

This  equation  is  called  a  "difference"  equation  because 
it  provides  an  update  to  the  previous  state  estimate  since  the 
observation  equations  are  non-linear.  If  they  were  linear 
there  would  be  no  need  to  iterate  and  the  optimal  least 
squares  estimate  is  just  ax  t  x,.  Nonetheless,  the  optimal 
estimate  for  the  non-linear  case  is  found  by  continually 
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updating  the  state  estinate  until  sone  convergence  criteria 
is  reached. 

Statistically, 


=  £(iX  SX^)  ,  (4.15) 

Where  E  is  the  expectation  operator.  Since  the  Taylor  series 
was  truncated  to  first  order, 

Sx  =  J  sy,  (4.16) 

where  j  is  the  Jacobian  of  x  with  respect  to  y.  Substitut¬ 
ing, 


=  £[Cr  Sy  (J  Syy] 

=  E(J  Sy  Sy^  CT).  (4.17) 

Since  J  is  a  constant  matrix,  it  can  be  brought  outside  of 
the  expectation  operator  producing 

=  J  £(  Sy  Sy^  ) 

=  cr  E>,  (4.18) 

Solving  for  IP^, 

=  a"‘  (cr")-‘.  (4.19) 

Since  J  is  an  orthogonal  symmetric  matrix,  s  J  \ 

(J).  (4.20) 
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Therefore,  given  the  functional  relationship  between  x  and  y, 
any  covariance  matrix  in  one  space  can  be  transformed  into 
another  space  via  the  above  association. 

VDobseryabil i ty 

One  of  the  main  problems  with  least  squares  is  that  the 
initial  estimate  must  be  sufficiently  "close''  to  a  solution 
to  ensure  convergence.  This  is  due  to  the  fact  that  the 
Taylor  series  was  expanded  about  a  reference  orbit  and  trun¬ 
cated  beyond  the  linear  term.  The  Gauss  orbit  determination 
method  produces  the  initial  state  estimate  for  input  into  the 
least  squares  estimator.  This  estimate  is  a  three  dimensional 
estimate,  however  the  data  being  fit  is  two-dimensional.  In 
other  words,  one  can  determine  the  direction  to  the  target 
with  respect  to  the  sensor,  but  the  data  contains  no  range 
information.  Thus  there  will  be  an  unobservable  eigenvector 
parallel  to  the  slant  range  vector  corresponding  to  an  infi¬ 
nite  eigenvalue  for  IE*. 

In  order  for  an  estimate  to  exist,  the  matrix  H^WH 
must  be  invertible.  This  is  called  the  observability  condi¬ 
tion  (19:58).  For  our  case,  the  zero  eigenvalue  of  pre¬ 
vents  inversion  of  this  matrix  to  form  the  covariance  but 
since  computers  introduce  roundoff  and  3ther  errors  we  can 
invert  this  ill-conditioned  matrix.  If  this  matrix  is  used 


41 


to  find  the  state  estimate  the  solution  will  rapidly  diverge 
due  to  an  uncontrollable  tx.  To  control  Ax,  decompose 
Into  its  eigenvalues  and  determine  the  contribution  of  the 
unobservable  eigenvector  Into  the  state  update. 

Let  T*  be  the  matrix  of  eigenvectors  of  1E*'\  so 

T  =  (  C,  }  (4.21) 

and  the  decomposed  P  '  matrix  Is  therefore 

P*‘=--TAT-‘.  (4.22) 

where  A  Is  the  diagonal  matrix  of  eigenvalues  x,  correspond¬ 
ing  to  ^1.  Taking  the  least  squares  equation  and  bringing  the 
covariance  matrix  to  the  left  hand  side 

P“  AX  »  H^H.  (4.23) 

Substituting  the  decomposed  matrix  and  premultiplying  each 
side  by  'Z*'^  produces 


A  T'‘AX  =•  (4.24) 

Letting  Ay  =  T'^AX  transforms  the  state  space  correction 
vector  Into  the  elgenspace  spanned  by  the  vectors  of  T'.  Upon 
examining  a  ,  there  will  be  a  zero  eigenvalue  and  the  corre¬ 
sponding  element  of  Ay  Is  its  contribution  to  the  correction 
vector  In  the  state  space.  Setting  this  element  of  Ay  to  zero 
eliminates  this  contribution ,  hence  transforming  this  vector 
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back  to  AX  produces  a  new  correction  vector  which  will  produce 
convergent  updates. 

Performing  this  procedure  is  not  without  its  drawbacks. 
Since  the  correction  contribution  along  the  slant  range  vector 
was  nulled,  there  are  many  different  solutions  to  the  azimuth 
and  elevation  data  fit.  Hence  any  reasonable  input  slant 
range  will  produce  random  residuals  upon  convergence  though 
the  slant  range  derived  from  the  converged  solution  will  not 
vary  much  from  the  input  value. 

Range  Determination 

Based  on  the  data  alone,  the  slant  range  from  the  sensor 
to  the  target  cannot  be  realized  since  the  information  con¬ 
tained  in  the  data  set  is  purely  direction  with  respect  to  the 
orbiting  platform.  The  best  range  estimate  one  can  obtain 
given  the  data  is  from  the  Gauss  method.  Even  then  the  ana¬ 
lyst  must  know  something  about  the  physics  of  the  problem  and 
the  geometry  involved.  However,  there  is  another  inherent 
piece  of  information  which  can  be  derived  from  the  data  —  the 
rate  at  which  the  target  passes  through  the  field  of  view. 

Assuming  the  platforms  and  objects  are  in  circular  co- 
planar  orbits,  the  angular  rate  at  which  satellites  revolve 
about  the  earth  is  the  mean  motion,  n,  in  radians  per  second. 
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r 


(4.25) 


n 


a’ 


i/» 


where  a  Is  the  senimajor  axis  of  the  orbit  and  (i  is  the 
Earth's  gravitational  parameter.  Defining  w,  as  the  angular 
rate  of  the  platform  and  a,  its  semimajor  axis,  the  target 
angular  rate  u  can  be  expressed  in  terms  of  a  variation  £a 
away  from  a,,  or 


u> 


'./■> 


I  (a,  +  Sa)^  J 


(4.26) 


The  above  equation  can  be  rewritten  in  a  binomial  series 
in  terms  of  the  ratio  ia/a. 


(j) 


00 


n«0 


(3/2)„(-l)" 

n! 


sa  y 
a,  . 


(4.27) 


where  the  notation  (a)^  detotes  the  Pochhammer  operator*. 
Assuming  6a  is  much  less  than  a,,  the  third  order  series  re¬ 
presentation  is 

w  «  w,  {  1  -  3/2a  +  15/8a’  -  35/16a’  )  ,  a  =  6a/a,.  (4.28) 

Therefore  the  angular  rate  of  the  target  relative  to  the 
sensor  is 


•The  Pochhammer  operator  (a)„  =  (a) (a+1) (a+2) . . . (a+n-l) . 
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(w  -  U),)  =  W,  {  -  3/2a  +  15/8o’  -  35/16a'  ), 

a  =  6a/a,.  (4.29) 

Since  both  objects  are  assumed  in  circular  orbits,  the 
radial  velocity  is  equal  to  zero.  This  makes  the  velocity 

merely  the  radius  times  the  angular  rate  in  radians  per  sec¬ 

ond,  or  v=ru.  Substituting  the  seroimajor  axis  for  radius 
makes  the  velocities  become 

v  =  u!  (a,  +  «a)  ,  V,  =  w  a,  (4.30) 

and  the  relative  inertial  target  velocity 

($v  =  (V  -  V,)  =  a,  (w  +  Lj,)  -  ia  w 

=  6a  w,  (  -  1/2  +  3/8a  -  5/160’  +  35/1280’  ), 

o  =  6a/a,.  (4.31) 

The  target  position  with  respect  to  the  platform  is  given 
by  the  vector  equation 


r  =  p  +  R, 

r  being  the  target,  R  the  platform  and  p  the  relative  position 
(slant  range)  vectors.  Since  p  is  measured  in  the  sensor 
coordinate  frame,  the  target  velocity  r'  becomes 

r'=  R'  +  p'  +  «  X  p  (4.32) 


45 


where  '  indicates  the  time  derivative.  Moving  to  the  left 
hand  side  and  expanding  p  into  its  components, 

r'  ”  *  iv  =  ‘*'®i  X  p®,.  (4.33) 

The  earth  subtends  approximately  17*  at  geosynchronous 
altitude.  A  simulation  of  two  coplanar  circular  orbits  and 
determing  the  relative  motion  of  a  lower  altitude  satellite 
with  respect  to  a  near-geosynchronous  satellite  has  shown  that 
within  the  viewing  cone,  the  sensor  relative  aspect  angle  rate 
9>'  is  constant.  Though  6'  is  constant,  the  value  of  fi'  varies 
with  slant  range  distance. 

Assuming  p  is  constant  through  the  field  of  view,  and 
since  6a  is  small,  p  »  6a,  the  relative  inertial  velocity  can 
be  written  in  terras  of  tangential  components  only,  or 

6v  =  6a  (  B'  +  u  ) .  (4.34) 

This  makes  the  expression 

B'  =  {  -3/2  +  3/8a  -  5/16o’  +  35/120a’  ;  , 

a  =  6a/a,  (4.36) 

which  can  be  solved  for  a  since  B'  is  known  from  the  data. 
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V.  Results 

Three  events  were  analyzed  with  the  developed  program, 
two  of  which  had  single  data  set  collections.  The  third  event 
was  quite  unique  in  that  there  were  four  consecutive  days 
worth  of  data  collected  on  the  same  object.  This  multi-day 
collection  was  verified  by  analyzing  the  intensity  profile 
history  for  each  collection.  Each  data  set  had  the  identical 
rotation  rate,  nearly  the  same  intensity  level,  and  started 
approximately  72  seconds  earlier  each  successive  day  pinpoint¬ 
ing  6a.  The  focal  plane  traces  are  shown  in  Figure  5.  For 
classification  purposes,  the  single  data  set  events  will  be 
designated  Event  A  and  Event  B.  The  multiple  event  will 
merely  be  called  the  multi-day  event. 

Event  A  and  B  Analysis  Results 

By  analyzing  Events  A  and  B,  deriving  the  slant  range 
from  the  angular  rate  is  plausible  only  when  the  ratio  ja/a, 
is  sufficiently  large.  Figure  6  illustrates  the  simulated 
relative  viewing  angle  from  a  spaceborne  platform  of  a  co- 
planar  target  in  a  circular  orbit.  For  ranges  close  to  the 
sensor  (within  1000  km),  6  has  constant  slope  -3/2w,  within 
the  viewing  cone  produced  by  the  limb  of  the  earth  at  geo¬ 
synchronous  altitude  (8.5"  or  .148353  radian  half-an¬ 
gle).  Figure  7  shows  that  for  small  6a/a,,  increases 

linearly  with  the  range  ratio  p/a,  at  a  rate  of  .375  to  l 

47 


I 


Figure  5 single  Data  Set  Focal  Plane  Traces 

which  was  expected  from  the  expansion  (realistically,  ^a/a,  < 
0.024,  which  is  an  optimistic  figure).  This  means  the  sum 
3/2w,  +  6'  roust  lie  between  -0.008877  and  zero  to  yield  a 
semi-major  axis  change  of  less  than  1000  km,  implying  that 
B'<0  since  the  change  in  semi-major  axis,  and  hence  slant 
range  must  be  negative  due  to  the  equation's  derivation  (ex¬ 
pansion  of  the  quantity  Sol  +  a,).  These  constraints  on  0' 
gave  rise  to  another  problem. 

The  constant  term  in  the  series  equals  1.09079624  X  10’* 
radians  per  second,  which  is  approaching  the  data  noise  level. 
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Figure  6  Simulated  View  Angle  Rates  for  Various  Targets 


Figure  7  Actual  Angular  Rate  Ratio  vs  Range  Ratio 
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Hence,  the  random  noise  on  5*  easily  corrupts  the  sum 
3/2  B*  producing  unreasonable  range  results  since  the  sum 
must  be  multiplied  by  8/3a,  (~  112645  km)  for  a  first  order 
approximation.  Therefore,  the  slant  range  calculation  is 
extremely  sensitive  to  data  corruption,  making  this  method  of 
range  determination  unfeasible.  This  result  is  in  agreement 
with  Hrastar  who  concluded  that  even  with  three  onboard  sen¬ 
sors,  small  uncertainties  in  target  location  produces  large 
uncertainties  in  the  range  (11:3). 

The  Gauss  orbit  determination  method  worked  fairly  well 
for  determining  an  initial  orbit  and  a  relationship  between 
r  and  v.  Some  of  the  derived  range  solutions  were  unrealistic 
since  they  produced  targets  either  behind  the  sensor  or  too 
far  away  to  be  detected.  Noise  was  a  minimal  contributor  to 
this  problem  since  the  data  were  smoothed  via  an  analyst 
supplied  n“  order  least  squares  polynomial.  In  fact,  chang¬ 
ing  the  polynomial  order  produced  negligible  change  in  the 
computed  radius  vector  through  3'“  order. 

After  fitting  the  data  on  these  events  with  an  arbitrary 
slant  range,  the  mean  of  the  residuals  was  on  the  same  order 
as  the  standard  deviation,  indicating  a  bias  existed.  This 
bias  essentially  shifted  the  estimated  state  vector  away  from 
the  nominal  solution  producing  unreasonable  target  selection 
choices . 
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Each  targets'  data  set,  regardless  of  its  size,  spans  at 
best  approximately  3%  of  its  orbit  on  a  single  pass  since  the 
target  is  also  near  geosynchronous  altitude.  The  slant  range, 
which  can  be  approximated  by  a  change  in  semi-major  axis 
between  the  target  and  sensor,  is  unknown;  only  the  direction 
from  the  sensor  at  any  one  time  is  known  given  the  data.  For 
both  of  the  single  set  events,  the  least  observable  eigen¬ 
vector  was  exactly  aligned  with  the  direction  cosine  vector 
from  the  sensor  to  the  target  with  a  near-infinite  eigenvalue. 
Therefore,  each  slant  range  chosen  by  the  analyst  will  fit  the 
data  equally  well.  Figure  8  through  Figure  11  show  the  resid¬ 
uals  for  both  Event  A  and  Event  B.  Table  II  lists  the  data 
fit  statistics. 

The  equinoctial  element  set  was  computed  from  the  state 
vector  to  avoid  the  singularities  associated  with  the  other 
element  sets.  The  state  vector  covariance  matrix  was  also 
transformed  to  the  equinoctial  covariance  matrix  to  determine 
the  uncertainties  of  each  element,  with  the  Jacobian  being 
computed  via  numerical  partials.  Unsurprisingly,  all  of  the 
elements  with  the  exception  of  the  semi-major  axis  were  known. 
The  single  data  set  equinoctial  element  covariance  matrices 
for  a  given  range  are  shown  in  the  appendix. 

Since  the  sensor  position  is  assumed  exactly  known,  the 
smaller  the  slant  range,  the  smaller  the  semi-major  axis 
variance,  a,’.  Figure  12  shows  a  higher  order  growth  in  the 
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Figure  9  Event  A  Elevation  Residuals 
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Event  B  Elevation  Residuals 
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Table  II 


Single  Data  Set  Event  Converged  statistics 


Mean 

SL 

RMS 

Event  A 

Az: 

-3.72130E-8 

1.38461E-4 

1.37496E-4 

El: 

-1.04324E-8 

2.77278E-5 

2.75345E-5 

Event  B 

Az: 

1.75770E-6 

5.25703E-5 

5.23875E-4 

El: 

-1.46506E-6 

1.67384E-5 

1.66801E-4 

semi-major  axis  certainty  the  further  the  target  is  from  the 
sensor.  This  makes  perfectly  good  sense  since  "£>,  = 
«(t,,tj)Pi  ♦(t,,ti)’.  In  fact,  from  this  relation,  o,'  should 
increase  in  the  order  of  p‘  at  best  since  ♦  is  at  most  second 
order . 

Multiple  Collection  Event 

It  was  quite  fortunate  that  a  series  of  collections  on 
the  same  object  was  located  in  the  data  base.  The  target  had 
a  periodic  intensity  profile  which  repeated  itself  upon  each 
track,  thus  ascertaining  that  it  was  the  same  object.  Each 
profile  consisted  of  a  series  of  specular  troughs  and  valleys 
indicating  that  it  was  multisurfaced.  A  rotation  rate  could 
not  be  derived  since  the  intensity  pattern  did  not  repeat  over 
the  duration  of  the  track.  Figure  13  shows  the  composite 
focal  plane  traces. 
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Knowing  the  change  In  track  start  times  from  the  previous 
day  would  give  the  analyst  an  idea  of  what  the  target  semima¬ 
jor  axis  should  be  based  on  the  sensor  period.  For  this 
event,  the  sensor  had  a  period  of  1436.154  minutes  and  each 
track  started  approximately  215  seconds  earlier  each  day  put¬ 
ting  the  target  70.5  km  away. 

since  two  body  motion  was  used  to  propagate  both  sensor 
and  target  state  vectors  over  three  days,  the  Integration  step 
size  had  to  be  decreased  a  significant  amount  (approximately 
60  second  Increments)  in  order  to  fit  the  position  ephemer- 
ides.  This  led  to  long  analysis  runs,  but  it  was  imperative 
that  each  succeeding  day's  sensor  epheroerides  (given)  match 
the  propagated  state  vector.  Also,  each  data  partial  had  to 
be  related  back  to  the  epoch  time  via  the  chain  rule. 

Using  the  a-priori  slant  range  value,  an  initial  orbit 
was  computed  from  the  Gauss  algorithm  for  the  first  day's 
worth  of  data.  A  trial  and  error  method  was  used  to  find  a 
slant  range  which  would  "hit”  the  next  day's  data.  Though  any 
range  would  fit  data  set,  only  one  slant  range  would  fit  two 
data  sets  collected  on  the  same  ob  ■  jt.  Each  latter  data  set 
became  easier  to  fit  once  the  first  two  data  sets  were  fit. 

Different  subroutines  were  "turned  off"  so  that  the  data 
could  drive  the  covariance  matrix  calculations.  For  instance, 
the  eigenvalue/observabJ lity  check  routine  was  not  needed 
since  the  covariance  matrix  was  now  invertible  —  the  slant 
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range  was  now  defined.  Data  partial  derivatives  were  not 
required  between  collections.  Hence  that  routine  was  bypas- 
sed  until  the  first  data  point  of  the  next  set  was  reached. 

The  entire  four  data  sets  were  fit  using  the  nominal 
Gaussian  standard  deviations  listed  in  the  sensor  specifica¬ 
tions.  As  expected,  the  semi-major  axis  error  dropped  signif¬ 
icantly  to  within  4  km.  The  error  for  the  remainder  of  the 
elements  also  was  reduced,  but  not  to  the  extent  as  the  semi¬ 
major  axis  error.  Figure  14  and  Figure  15  show  the  multiple 
day  event  data  fit  residuals.  Note  the  end  of  each  data  set 
after  careful  inspection  of  the  plots. 

Table  III  lists  the  derived  classical  element  set  and  the 
respective  element  errors.  Any  element  which  involves  veloci¬ 
ty  in  its  computation  (i.e.  eccentricity,  argument  of  perigee, 
mean  anomaly  and  semi-major  axis)  has  a  higher  uncertainty 
than  the  pure  position  derived  elements  (ascending  node  and 
inclination)  since  error  is  extremely  sensicive  to  changes  in 
position  derivatives. 

Once  the  data  was  fit,  the  obtained  fastwalker  state 
vector  was  propagated  to  find  the  range  at  the  nadir  crossing 
for  each  of  the  four  days,  shown  in  Table  IV.  The  target 
appeared  to  have  drifted  approximately  40  km  further  from  the 
sensor  each  day  for  two  days  and  then  stopped.  This  effect 
was  due  to  the  relative  changes  in  the  orbit  between  the 
sensor  and  target  since  the  eccentricity  is  not  exactly  zero. 
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Figure  14  Multiple  Day  Event  Azimuth  Residuals 


Figure  15  Multiple  Day  Event  Elevation  Residuals 


Tabid  III  Derived  Multiple  Day  Event  Results 


Semi-major  Axis,  a 

42167.15  km  a  = 

3.6  Km 

Eccentricity,  e 

.00298  a  ~ 

.000213 

Inclination,  i 

5.664“  a  - 

.000305“ 

Ascending  Node,  n 

70.566“  a  = 

.000746“ 

Argument  of  Perigee,  w 

316.800“  a  = 

.0198“ 

Mean  Anomaly,  M 

96.489“  a  = 

.00281“ 

Azimuth 

Mean:  8.17632E-7 

RMS:  4.96381E-5  a:  4 

.99205E-5 

Elevation 

Mean:  5.61084E-7 

RMS:  1.17450E-5  a:  1 

.20884E-5 

Table  IV  Multiple  Day  Event  Nadir  Crossing  Ranges 


Day  0; 

79.0  km 

Day  2 

176.2 

km 

Day  !■ 

130  km 

Day  3 

180.4 

km 
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Fastwalker  orbit  determination  cannot  be  performed  with 
a  single  data  set.  Multiple  simultaneous  collections  by 
different  sensors  or  another  sighting  by  the  same  sensor  on 
a  different  day  is  required  to  determine  the  orbital  element 
set  with  any  certainty .  Hence  any  future  occurrences  with 
other  censors  cannot  be  computed  accurately  since  errors  are 
also  jpagated  along  with  the  state. 

A  near-perfect  estimate  of  the  sensor  position  and  veloc¬ 
ity  is  required  so  any  errors  contributed  to  the  results  will 
be  due  t''  the  data  and  not  the  ephemeris.  This  is  accom¬ 
plished  by  assuming  the  first  ephemeris  position  point  is 
known  and  fitting  the  remainder  of  the  position  ophemerides 
by  estimating  the  velocity.  This  adds  confidence  to  the 
propagated  state  vector  when  fitting  multiple  day  colJections. 

The  Gauss  method  provides  a  good  initial  orbit  provided 
the  data  does  not  have  any  biases  (zero  mean),  especially 
since  the  method  assumes  perfect  data.  Hence,  no  error  sta¬ 
tistics  are  associated  with  the  elements.  Once  the  eighth- 
order  polynomial  is  factored,  the  method  occasionally  produces 
possible  target  ranges  where  the  target  ^^ould  not  be  detected 
due  to  the  lower  reflected  sunlight  intensity  level  or  a 
negative  range  meaning  the  target  is  behind  the  sensor.  These 
outputs  are  naturally  impossible  and  the  method  fails.  The 
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method  does  however  provide  a  relationship  between  the  target 
position  and  velocity  vectors  as  a  function  of  range. 

Since  the  target  is  essentially  coplanar  with  the  sensor, 
the  range  is  the  difference  in  the  target  semi-major  axis  from 
the  sensor.  Unfortunately,  the  least  observable  eigenvector 
derived  from  the  covariance  matrix  is  along  the  line  of  sight 
indicating  that  any  range  will  fit  the  azimuth  and  elevation 
data  equally  well.  The  least  confident  element  is  the  semi¬ 
major  axis;  its  error  increases  as  a  function  of  p*.  The 
balance  of  the  element  errors,  whether  they  be  classical  or 
equinoctial,  increase  at  a  much  slower  rate.  For  data  sets 
with  T\f  '  jro  mean,  the  remainder  of  the  elements  are  confi¬ 
dently  known  when  the  target  is  within  100  km.  Otherwise  the 
errors  grow  almost  as  fast  as  the  semi-major  axis  error. 

Multiple  sightings  on  different  days  require  a  relatively 
small  integration  step  if  two  body  motion  is  to  be  used  for 
propagating  the  sensor  and  initial  fastwalker  state  vectors. 
The  serai-iuajor  axis  error  is  significantly  reduced  with  each 
additional  data  set.  Once  two  data  sets  are  fit  (which  is  a 
chore  in  and  of  itself),  the  third  and  later  sets  become 
easier  to  analyze  since  the  epoch  estimate  is  more  refined. 
The  major  problem  is  finding  a  range  which  will  produce  an 
orbit  consistent  with  the  second  data  set.  Perturbations  do 
not  play  a  major  role  since  these  effects  are  insignificant 
over  the  time  period  of  the  analyzed  events  and  they  are  also 


61 


absorbed  in  the  data  noise.  Though  the  data  residuals  still 
have  trends,  the  RMS  values  are  acceptible  given  the  bias 
problems . 

Care  must  be  taken  to  ensure  that  the  same  object  is 
being  tracked  else  many  hours  will  be  wasted.  This  is  done 
by  comparing  intensity  history  patterns  of  various  detected 
objects.  The  intensity  level  is  not  important  since  it  is  a 
function  of  the  materials  used  to  build  the  target,  location 
of  the  sun  at  the  time  of  collection,  target  attitude  with 
respect  to  the  sensor  and  other  effects  in  addition  to  dis¬ 
tance  from  the  sensor.  Only  the  specular  pattern  is  important 
since  most  of  the  satellites  in  geosynchronous  orbit  are 
spinners  with  despun  shelves.  This  pattern  gives  the  analyst 
an  idea  of  the  target  shape  and  the  spin  rate.  It  is  similar 
to  a  fingerprint  —  an  identifying  trait  which  discerns  one 
object  from  another. 
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VII.  Regojaroendatiflna 


The  topic  of  this  thesis  has  not  been  deeply  investigated 
by  many  analysts  within  the  space  and  intelligence  community. 
It  is  only  recently  that  fastwalker  analysis  has  become  an 
important  issue  due  to  an  ever  increasing  number  of  geosyn¬ 
chronous  satellites  and  the  inability  to  track  and  catalog 
these  objects  with  ground-based  sensors.  Therefore  many 
related  topics  are  available  for  research. 

The  main  hurdle  to  overcome  with  fastwalker  analysis  with 
a  single  data  set  is  the  range  ambiguity.  Rederiving  the 
Gauss  orbit  determination  method  to  include  state  vector  error 
statistics  will  be  able  to  provide  the  estimator  information 
about  what  it  is  calculating.  Since  the  direction  cosine 
vectors  are  fit  with  a  least  squares  polynomial;  there  is  some 
error  information  present  which  could  be  related  to  the  ini¬ 
tial  state  estimate  from  the  Gauss  algorithm.  It  is  possible 
to  transform  these  error  statistics  into  some  type  of  a-priori 
covariance  matrix  which  in  turn  can  be  added  to  the  least 
squares  covariance  matrix.  This  will  control  the  estimation 
process  iteration  and  reduce  the  need  to  decompose  the  least 
squares  covariance  matrix. 

Secondly;  devising  a  method  to  remove  biases  from  the 
data  will  help  allow  the  Gauss  method  to  perform  better  with 
the  given  data  thus  providing  a  more  confident  initial  state 
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vector  to  the  estimator.  This  can  be  combined  with  the  above 
recommendation  to  again  furnish  the  estimator  with  error 
statistics  about  what  it  is  receiving. 

Lastly,  one  can  perform  intensity  history  profile  compar¬ 
isons  to  determine  when  the  same  object  is  detected  by  two 
orbiting  sensors  in  different  locations.  This  will  involve 
perturbation  theory  and  propagating  state  vectors  (or  orbital 
elements)  over  long  time  periods  and  then  fit  two  or  more  data 
sets.  Accurate  propagation  models  using  the  HAMING  integrator 
will  be  required  (19:B-12). 


64 


Appendix:  Covariance  Matrices  for  Events  A  and  B 


The  following  matrices  are  typical  single  data  set  event 
covariance  matrices.  Each  matrix  is  for  a  specific  event  and 
the  matrix  elements  are  functions  of  range  and  how  well  the 
data  was  fit.  The  given  ranges  were  derived  by  the  Gauss 
algorithm.  Read  the  elements  from  left  to  right  and  top  to 
bottom  as  follows:  a,  h,  k,  Xg,  p,  q  (i.e.  is  the  partial 
derivative  of  h  with  respect  to  k). 

Event  A  Equinoctial  Element  Covariance  Matrix  at  t..  o  =  25km 

0.46298E+00  0.36882E-03  0.60340E-03  0.41122E-04  0.12580E-02  0.11407E-02 
0.36882E-03  0.29422E*06  0.48L27E-06  0.32800E-07  0.10034E-05  0.90989E-06 
0.60340E-03  0.48127E-06  0.78739E-06  0.53660E-07  0.16415E-05  0.14885E-05 
0.41122E-04  0.32800E-07  0.53660E-07  0.36570E-08  0.11187E-06  0.10144E-06 
0.12580E-02  0.10034E-05  0.16415E-05  0.11187E-06  0.34222E-05  0.31032E-05 
0.11407E-02  0.90989E-06  0.14885E-05  0.10144E-06  0.31032E-05  0.28140E-03 

Event  B  Equinoctial  Element  Covariance  Matrix  at  t„.  p  =  425kro 

0.20772E+07  0.87512E+01  0.58293E+01  0.32037E+01  0 . 21847E+02 -0 . 56257E+02 
0.87512E+01  0.37038E-04  0.24631E-04  0.13496E-04  0 . 92037E-04-0. 23700E-03 
0.58293E+01  0.24631E-04  0.16394E-04  0.89899E-05  0 . 61306E-04-0 . 15787E-03 
0.32037E+01  0.13496E-04  0.89899E-05  0.49410E-05  0 . 33694E-04-0 . 86766E-04 
0.21847E+02  0.92037E-04  0.61306E-04  0.33694E-04  0 . 229 77E-03-0 . 59168E-03 
-0.56257E+02-0.23700E-03-0, 15787E-03-0.86766E-04-0.59168E-03  0.15236E-02 
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